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We present a minimalistic approach to simulations of force 
transmission through granular systems. We start from a con- 
figuration containing cohesive (tensile) contact forces and use 
an adaptive procedure to find the stable configuration with no 
tensile contact forces. The procedure works by sequentially 
removing and adding individual contacts between adjacent 
beads, while the bead positions are not modified. In a series 
of two-dimensional realizations, the resulting force networks 
are shown to satisfy a linear constraint among the three com- 
ponents of average stress, as anticipated by recent theories. 
The coefficients in the linear constraint remain nearly con- 
stant for a range of shear loadings up to about .6 of the nor- 
mal loading. The spatial distribution of contact forces shows 
strong concentration along "force chains". The probability 
of contact forces of magnitude / shows an exponential falloff 
with /. The response to a local perturbing force is concen- 
trated along two characteristic rays directed downward and 
laterally. 

PACS numbers: 46.10.+Z, 83.70.Fn 



I. INTRODUCTION 

The fragility of granular matter is a longstanding pre- 
occupation of engineers ||l|,^ and a recent preoccupation 
of physicists pB]. By granular matter we mean a static 
assembly of hard, spheroidal grains whose contact forces 
may be compressive but not tensile. Thus granular mat- 
ter is noncohesive. Mohr and Coulomb recognized a 
fundamental continuum consequence of the noncohesive 
state. There can be no co-ordinate system in which the 
shear stress exceeds some fixed multiple fi of the normal 
stress M. This "Mohr-Coulomb" condition limits the 
stresses that a granular material can support, and thus 
amounts to a form of fragility. When this condition is 
violated, building foundations settle and embankments 
slip. 

Modern civil engineering practice g views the stress 
field in a granular medium as divided into elastic and 
plastic zones. The stresses in these plastic zones are at 
the Mohr-Coulomb limit, and thus these zones are at 
the margin of stability. The stress in the elastic zones is 
within the bounds of stability and thus the stress here is 
transmitted as in an elastic body. 

Recently attention has turned to the microscopic ori- 
gin of the macroscopic fragility of granular media. The 
microscopic pattern of contact forces and bead motions 
shows strong local heterogeneity and history dependence 
|^-|ll|. The history of prior motion in a region clearly 



influences the way it transmits forces. The prior motion 
may affect the /i coefficient in the Mohr-Coulomb law, 
the elasticity tensor, or further constitutive properties 
0. The question is, for a given history of relaxation to a 
static state, how are these forces transmitted and what 
range of forces can be supported. The transmission of 
forces can be expressed as a linear-response property of 
a granular pack. An infinitessimal force is added to the 
bead at position xq and the corresponding incremental 
force on a contact at r is determined. For sufficiently 
small perturbations of a finite pack, this linear response 
function G is well defined. It depends on the shapes and 
sizes of the beads, their frictional properties, and how the 
pack was constructed. 

If the perturbing forces become too large, motion oc- 
curs. Beads shift their positions and form new contacts. 
This motion may be reversible, so that the beads re- 
turn to their original positions when the perturbation 
is removed. This motion may also be irreversible, with 
the positions altered after removal of the perturbation. 
The thresholds for reversible and for irreversible motion 
are fundamental ways to characterize the nonlinear re- 
sponse of the pack. For any given pack there is a weakest 
perturbing force distribution that causes motion. This 
threshold force may go to zero as the size of the pack 
grows. Many simulations have sought to characterize 
the above features of a granular pack [p2|-14|. These 
studies model the system in a realistic way that requires 
detailed specifications and many parameters. This de- 
tail makes it difficult to discern which observed features 
are inescapable consequences of the granular state, and 
which are properties of the particular realization. In this 
study we take the opposite approach, sacrificing realism 
for the sake of simplicity. We seek the simplest system 
that shows the instabilities of noncohesive material. Thus 
our system consists of frictionless, spherical beads, which 
have been deposited into a container one at a time and 
not moved thereafter. Such a system develops tensile 
contacts. 

To avoid these contacts, we must define some motion 
that evolves the pack to a more stable state. Again we 
choose a procedure favoring simplicity rather than real- 
ism. We seek the stable state attainable with minimal 
disturbance from the initial state. Accordingly our pro- 
cedure does not move the beads, but rather removes and 
adds contacts one at a time in order to attain a stable 
contact network. In this sense our simulation is an adap- 
tive network. 

This network is isostatic pq]: the contact forces are 



determined from the applied forces solely through the 
force equilibrium of each bead, without reference to bead 
displacements or material deformation. 

Our adaptive method demonstrates that frictionless 
granular materials can be mechanically robust. For a 
given load, the simulation converges to a state of no 
tensile contacts. A change in the applied load of order 
unity can be applied with only minor shifts in the con- 
tacts. Further, the three components of the stress in two 
dimensions obey a constitutive law of the "null stress" 
type: a weighted sum of the three components vanishes, 
the weights depending on the packing but not on the 
loading. 



II. RESPONSE FUNCTION 

Our system is a set of spherical beads, whose radii are 
chosen randomly within a moderate range. These beads 
are supported on one side, called the bottom, with a layer 
of fixed spheres. The width w of the system is much 
larger than a bead. The beads are arranged densely in 
this space up to a height h. A fixed downward force Fq 
is applied to each bead lying at the upper surface. We 
choose a configuration of beads that is mechanically sta- 
ble: the normal forces acting at the bead contacts oppose 
the applied forces Fq and prevent motion. Initially, we 
allow for tensile (negative) contact forces. 

We label the TV beads by an index a. Then we may 
denote the contact force from bead a to bead (3 by the 
scalar fa/i- Newton's Third Law dictates that for all a, 
P, fj3a = fap- The Nc contact forces are constrained 
by the requirement that the total vector force on each 
bead vanish. In d dimensions, there are evidently d N 
such constraints. The bead positions Xq, are likewise con- 
strained by the geometrical condition that the distance 
between two contacting beads a and (3 must be the sum 
of their radii r^ -\- rp. There are N^ such constraints for 
the d N quantities Xq. If all these constraints are inde- 
pendent, the number N^ of contacts must be exactly d N 
p6[ . Then the system is isostatic: the d N force balance 
equations are just sufficient to determine the Nc contact 
forces. The equations of force balance may be written 



/ ^ fapflafj — Fq 
/3(a) 



(1) 



Here (3{a) denotes the set of contacting neighbors of bead 
a, Fq, is an external force applied to this bead, fap and 
the unit vector hap represent the magnitude and (fixed) 
direction of the contact force between the two beads a 
and (3. Since all the above equations are linear, the re- 
sponse of the system to a given external forcing is deter- 
mined by the response function G: 

lap = G(a/3|7) • F^ (2) 

The response function G determines not only the re- 
sponse to an external force but also the global displace- 
ment field associated with local geometrical perturbation 



of the network. In order to see this, let us assume that 
the packing is subjected to external forcing F^. Then 
we relax exactly one of the N geometric constraints, and 
change infinitesimally the distance between two contact- 
ing beads, r^f} = |x^ — Xq|. As long as the connectiv- 
ity of the network does not change, its motion is non- 
dissipative. This means that the work done to distort 
the packing locally, Srapfap is the work against external 
forces, i. e. 

(5x^ • F-y = Srapfap = fra/3G(a/3|7) • F-^ (3) 

The above equation should be valid for any set of ex- 
ternal forces; hence. 



6xj = G{a(3W)6rai3 



(4) 



We conclude that G is the response function both for 
contact force and the displacement field. Note that the 
displacement discussed here is not due to deformation of 
the beads. It corresponds to a "soft mode" that pre- 
serves the distances between all contacting beads other 
than the perturbed contact. 

In a general case, finding the response function for a 
given configuration is a non-local problem, which requires 
solving a set of linear Eq. (|l|) . The task becomes much 
easier for the case of sequential packing. This is created 
by adding one bead at a time. The requirement of me- 
chanical stability implies that any newly-added bead has 
exactly d "supporting" contacts (in d-dimensional space). 
If all the contacts were permanent and this d-branch tree 
structure were not perturbed by the future manipula- 
tions, the response function might be found by a simple 
unidirectional projection procedure. Indeed, since there 
are exactly d supporting contacts for any bead in a se- 
quential packing, the total force F^ including external 
force Fq and that applied from the supported beads can 
be uniquely decomposed onto the corresponding d com- 
ponents, directed along the supporting unit vectors n^^. 
This gives the values of the supporting forces. The /'s 
may be compactly expressed in terms of a generalized 
scalar product (...|...)^: 



fap 



T-aP 



(5) 



The scalar product (...|...)^ is defined such that 
{hap\hapi) a = ^(30'- for the supporting contacts /3, /3' 
of bead a. In general, it does not coincide with the con- 
ventional scalar product. The resulting response func- 
tion, G{a(3\j), can be calculated as the superposition of 
all the projection sequences (i.e. trajectories), which lead 
from bead 7 to the bond a/3: 



G(a/3|7) = 2^ \h^a^)^{n'■fa^\na,a2)a,■■■{'^aka\nap)a 

(7^Ql...— >0!^/3) 

(6) 

Here the summation is done over all the trajectories 
(7 — > Qfi — > ... ^ ak -^ a ^ (3) such that any bead 
in the sequence is a supporting neighbor of the previous 



III. ADAPTIVE NETWORK SIMULATION. 

For a large enough system, sequential packing is not 
compatible with the requirement of non-tensile contacts. 
Anytime when this requirement is violated, a rearrange- 
ment occurs and system finds a "better" configuration. 
One might expect that this would make the problem of 
force propagation a dynamic one. However, it is possi- 
ble to limit oneself to a purely geometrical consideration, 
following the ideas of the previous section. 

Suppose afi is a "bad bond", whose contact force is 
negative (tensile). This means that the network would 
move in such a way that the two beads, a and (3 are 
taken apart. In other words, the soft mode associated 
with the perturbation of a(3 bond is activated, and for 
small enough displacements all the beads move in accor- 
dance with Eq. (||). The motion stops when a replace- 
ment contact is created, i.e. when a gap between any two 
neighboring beads closes. 

In this work, we limit ourselves to this linear approx- 
imation. It should be understood that Eq. (H) is cor- 
rect only for infinitesimal displacements, and in a general 
case one should account for the evolution of the response 
function in the course of the rearrangement. We avoid 
the problem of changing G by permitting only infinites- 
simal motion in the model. We imagine that the "bad 
bond" gets deactivated, and it is replaced with a rigid 
"strut" between two neighboring beads that were not in 
contact in the previous configuration. There is a natu- 
ral choice for where the strut should be placed to cause 
minimal disturbance. Each pair of non-contacting neigh- 
bors 7(5 has a gap r^s ~ T-y — fs- When the contact a/? 
is removed, the distance rap is allowed to change; this 
change alters the gaps of other neighbors ^5 as specified 
by Eq. |4[ Extrapolating this linear-response equation, 
motion Srap required to close gap ^5 is given by 



5ra 



P 



r^s 



rs 



fi^s- [G(a/3 7) - G(a/3 5)] 



(7) 



(For many choices of ^5 the required 5rap is infinite since 
the a(i contact has no effect on the ^5 gap.) Using this 
formula, we identify the gap a' (3' which would require 
the smallest change of Tq,^ in order to close, and we link 
this pair by a strut. 

After we have found the replacement bond, the mod- 
ified response function can be found without solving the 
whole set of force balance equations (|l|) ! We denote the 
response function for the initial packing as Go- The new 
response function G must be such that there is no longer 
a contact force fap- In general there is such a force in 
the initial packing. However, we may alter this unwanted 
force by adding an external force to some other bead. 
We choose to add external forces to beads a' and j3' that 
mimic a contact force: the two forces are equal, opposite 
and directed along the unit vector nct'/3' joining them, 
with a strength denoted fa'p'- Our choice of the replace- 
ment pair guarantees that the effect of a' or /3' on the afi 



contact is non-vanishing. Then the force on this contact 
is given by 



/„^ = ;^F^.Go(a/3| 



7j 



(8) 



fa'li'[Go{al3\a') - Go(a/9|/3')] ■ fia'p' 



We may make this fap vanish by a proper choice of 
the external force, fa' (5' — '^^'P-y ■ G(a'/3'|7), where 
G(a'/9' 7) is the new response function, found by requir- 
ing that fap vanish in Eq. (||): 



G{a'l3'h) 



Go(a/3|7) 



fia'p' ■ [Go(a/3 a') - Go(a/3 /?')] 



(9) 



A contact force on an arbitrary contact is now deter- 
mined by a combination of external forces F^ and the 
above-determined fa'p'- This results in the following ex- 
pression for G{Xii\'-f) (A/i other than a'(3'); 



G(Am|7) = GoiXfi\j) 
-Goiaph) 



(10) 



fiq'p, ■ [Go{X^l\a') - Go(A/i|/30] 
n,,^, •[Go(a/3|a')-Go(a/3|/3')]' 



This prescription gives the response of the pack to a 
wide class of contact replacements. The prescription does 
not require either the initial or the final state to be stable: 
it allows tensile contacts. Using the contact replacement 
procedure we may investigate the stability of a pack sys- 
tematically. Our algorithm has two major stages: net- 
work preparation and its "mutation" via the contact re- 
placement scheme. By repeating this adaptive procedure 
sufficiently many times, one may hope to get the stable 
configuration without tensile forces (for a given loading) , 
just like the real system would do. There is a possibil- 
ity that the present geometry-preserving algorithm could 
not stabilize the network. For instance if the tangen- 
tial component of the surface force is strong enough, it 
is expected to initiate a macroscopic avalanche, as in a 
sandpile with slope exceeding the critical angle. This 
class of rearrangements is beyond the capabilities of our 
connectivity mutation scheme. This circumstance has 
even certain advantages: we can determine the critical 
slope from our simulations as the direction of the surface 
force at which the algorithm stops working. It should 
be emphasized that our algorithm can easily be mod- 
ified to incorporate the change in the network geome- 
try. The major reason why we use the above geometry- 
preserving ("strut") approximation is its much higher 
computational effectiveness. 



IV. SIMULATION DETAILS AND RESULTS 

A. Method 

We begin by creating a two-dimensional sequential 
pack of variable-sized discs by adding them one by one. 



The studied system has the following parameters: poly- 
dispersity 10% (bead radii from 1 to 1.1), number of 
beads, N from 250 to 500 (the major limitation is the 
computation time). Although there is no gravitational 
force acting on the beads in these simulations, the statis- 
tics of the packing can be varied by changing the "pseudo- 
gravity" direction, g. Namely, while adding a bead to 
the packing we require g to be directed between the two 
supporting contacts. Simultaneously, we calculate the re- 
sponse function, by using the sum-over-trajectories for- 
mula, Eq. (H). 

Then we apply certain load to beads on the surface. In 
the studied cases, the forces applied to the surface beads 
were all the same, with the only principal variable being 
the ratio of the two components fx/ fy [y is the vertical 
direction, the surface in average is parallel to x due to 
the periodic boundary conditions in the horizontal direc- 
tion). As long as the response function is given, we may 
find all the contact forces for a given load. As we have 
found, tensile contacts appear within a few beads from 
the surface. We analyze the sign of the contact forces one 
by one, from top to bottom (in the order opposite to the 
one in which the beads were originally deposited) . When 
a tensile contact is encountered, we follow the contact 
replacement procedure described in the previous section: 
find the new bond and modify the response function to 
account for the connectivity change. Now we repeat the 
procedure again starting from very top until there is no 
tensile contact left in the system Figure |l| shows a se- 
quence of four typical steps in this procedure. Evidently 
the removed and the added contact may be far apart. 




> 




FIG. 1. Four typical steps in our annealing sequence. The 
beads are pictured as circles and the bonds are shown as gray 
lines. Upper left to lower right shows moves 804-807. The 
entire annealing process required 1491 such moves. In each 
frame the removed contact is shown as a thick white line; the 
added contact is shown as a thick black line. 

It sometimes happens that this prescription does not 
remove the tensile contacts: the removal of a tensile con- 
tact continues to generate others. In this case we may 
modify our procedure for selecting the next tensile con- 
tact to remove. For example we may select the strongest 
tensile contact instead of that tensile contact having the 
largest sequence number. Such alternative prescriptions 
seem to have little effect on the force network, as dis- 
cussed in the next section. 



B. Variability and reproducibility 

While our bond replacement procedure mimics the way 
the real system should rearrange, our choice of the "bad 
bond" to be replaced is far more arbitrary. For instance, 
instead of checking the sign of the contact forces one by 
one from top to bottom, we could go the other way, or 
try to replace the contact bearing the largest negative 
force first. Neither of these prescriptions is very realis- 
tic; however, we find that the results are insensitive to 
the procedure. In order to probe this sensitivity, we 
compared the results from two different "annealing" pro- 
cedures. The first was the procedure described in the 
previous subsection. The second procedure is as follows. 
We perform exactly the same one-by-one check, as in the 
previous case, but do not remove a negative bond un- 
less the magnitude of the force exceeds certain tolerance 
threshold. When no more bonds remain in the system 
for a given threshold level, we reduce the tolerance and 
repeat the procedure. The threshold plays a role of tem- 
perature: keeping it finite allows us to deviate from the 
target (non-tensile) state of the system and explore its 
vicinity at the configuration space. The second an- 
nealing algorithm converges considerably faster than the 
zero-tolerance one. For instance, it took us from 1500 to 
3000 iterations to complete the original algorithm with 
500 beads, while the annealing procedure reduced the 
needed time to approximately 500-1000 steps. Interest- 
ingly, the variation of the convergence time is of order 
of the time itself. 

We have found that the contact configuration result- 
ing from the annealing procedure does differ from the one 
generated by the zero-tolerance algorithm (see Figure ||). 
However, we did not detect any statistically-significant 
variation of the ensemble-averaged properties of the final 
state obtained with the two methods. These properties 
included the average stress and the contact force proba- 
bility distribution function (PDF), presented below. 




FIG. 2. Pattern of contact forces from two different an- 
nealing runs on the same bead pack. Each contact force is 
shown as a hne joining the centers of the two contacting beads. 
Thickness of the hne is proportional to the magnitude of the 
force. Dark lines show forces obtained by replacing each ten- 
sile contact as it is encountered in a search from the top down. 
Light lines show the forces obtained when only strongly tensile 
contacts were replaced at first. Then when no further strongly 
tensile contacts remained, the rest of the tensile contacts were 
relaxed. The heavy horizontal lines along the bottom are ar- 
tifacts introduced by our rendering program. 



C. Macroscopic constitutive equation 



major hypothesis used in our earlier work Pq | to derive 
the constitutive equation of frictionless granular packing. 
One more hypothesis, used for derivation of the macro- 
scopic equation for stress is the mean field decoupling 
Ansatz. The average stress in a region of a sequential 
packing can be written ||lq| 






(11) 



^) 5(x„ - x) ( F„|n„^^ </3</3^"/3 ) : 



The sum /3 is over the beads that support the bead a. 
Our mean-field hypothesis consists in assuming that the 
force-related part of this average is independent of the 
geometrical part. We define 



f(x)^(^<5(x„-x)(Fj 



and 




</3<,3'^a/3 



(12) 



(13) 



Then our mean-field assumption amounts to the state- 
ment 

cr'^ = f • f (14) 

The mean field hypothesis can be directly checked for 
unperturbed sequential packing, were both fields f and 
f are well-defined. The results of such a check are rep- 
resented on Figure 0. The exact values of various stress 
components are shown to be in an excellent agreement 
with their evaluation based on the mean-field Ansatz. We 
conclude that the mean field is a very good approxima- 
tion at least for non-adaptive sequential packing. 



One of the crucial results of the simulation is that our 
geometry-preserving adaptive network algorithm does 
converge for a considerable range of force direction. It 
stops working when \fx/fy\ approaches 0.6 (for pack- 
ing prepared at vertical pseudogravity g. This suggests 
that the critical slope for the frictionless packing is about 
30 degrees, consistent with simple theoretical arguments 
and some experiments [|7| Note that this slope may con- 
siderably exceed the angle of repose in dynamic exper- 
iments and simulations because of hysteresis associated 
with the lack of damping in the frictionless system. Pre- 
sumably, the critical slope can be observed by quasi-static 
tilting of a zero-slope packing. 

Another interesting observation is that the eventual 
connectivity of the packing is not too different from the 
original sequential packing. For example, the 500-bead 
system needs up to 3000 iterative steps (rearrangements) 
to find the stable state, and yet only 150 out of 1000 con- 
tacts (15%) in the final configuration are different from 
the original network. This provides us with a solid back- 
ground for using the sequential packing as the zero-order 
approximation of the real network. This was one of the 



0.5 



^J^yy 




-0.5 ^ 
FIG. 3. Various components of stress tensor in original se- 
quential packing (before the adaptive stage), as functions of 
the direction of the applied force. Note a remarkable agree- 
ment between the mean field results (dashed lines) and the 
simulation data (solid lines). 
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4. The only free parameter of the stress tensor, 
yy, as a functions of the direction of the applied force, 



FIG 
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after the adaptive algorithm is completed. The solid lines 
show the simulations results corresponding to two different 
directions of "pseudo-gravity" vector, g. The dashed lines 
represent the null-stress law corresponding to the material 
tensor f computed for the original sequential packing. 

As long as the rearrangements are switched on, there is 
no obvious way to define the concept of supporting neigh- 
bor, and therefore the "force from subsequent beads", f 
is ill-defined as well. However, a more general meaning 
of constitutive Eq. ^Wj , is that the stress is parameteriz- 
able with some vector f , and the third-rank material ten- 
sor f establishes this parameterization. We now take f 
corresponding to the original pre-rearrangenient sequen- 
tial packing and probe the Ansatz by checking whether 
the total stress (after the adaptive procedure) can be ex- 
pressed as f ■ f . In other words, we compare the only un- 
known component of the stress a^x (two other are given 
by the boundary conditions) with its theoretical value ob- 
tained from f. We performed this check for two different 
classes of packing, corresponding to different directions 
of "pseudogravity" , gx/gy = and 0.2. The agreement 
between the two curves is surprisingly good as long as 
the direction of the applied force does not deviate too 
much from the preparation conditions (i. e. from the 
pseudogravity vector) , see Figure ^. 



a point source according to a wave- like equation |^ . In a 
medium where all non- vertical directions are equivalent, 
the force should propagate downward along slanting char- 
acteristic lines, whose slope is dictated by the only pa- 
rameter in the null-stress law. The responding region lies 
within the "light cone" bounded by these lines. In two 
dimensions, the response consists of two delta- functions 
traveling along the light cone. Disorder is expected to 
scatter the wave solutions of the pure system, thus re- 
sulting in a widening of the delta-peaks. This scattering 
could be sufficient to create qualitative new mesoscopic 
behavior from localization effects [Q. Our simulated 
system showed strong influences from disorder, as illus- 
trated in Figure o. Because there can be no vertical-force 
response at the top of the system, we observe a global 
anisotropy, with stronger responses below the source than 
above it. The response is also strongly heterogeneous. 

Our simulations allow us to perform ensemble averag- 
ing of the response stress field. Figure g shows the results 
of such averaging over 600 realizations of the network. As 
the perturbation propagates deep into the sample, the re- 
sponse function gets a two-peak shape, in a good agree- 
ment with the null-stress law. As expected, the peaks are 
broadened by the disorder, and one cannot resolve them 
immediately below the source. Another important ob- 
servation, which also supports the null-stress approach, 
is that the average response is virtually zero above the 
source. Note that we have studied only linear response 
of the system, so that the perturbation did not change 
the contact network. This need not to be the case in 
the experiments involving strong local perturbations M] 




FIG. 5. Two patterns of contact forces resulting from a 
perturbing point force in a single relaxed configuration. Per- 
turbation was an infinitessimal downward force applied at the 
points shown as heavy dots. A compressive contact force is 
indicated as a dark line joining the centers of the two contact- 
ing beads. Thickness is proportional to the magnitude of the 
force. Tensile forces are indicated as lighter gray lines. 



D. Response function 



As noted above, our system transmits forces in ac- 
cord with a null-stress constitutive property. Given the 
null-stress law, one may infer the corresponding response 
function G. The force transmission is transmitted from 



0.1 



/^ 



/ 

/ 
/ 


\ 

\ 
\ 


y=-9 

y=-3 


/ 

/ 

/0.06 - 


\ 

\ 
\ 
\ 


y=6 



yy 




-0.02 -L 

FIG. 6. ayy component of the ensemble-averaged response 
to a unit vertical force applied at the origin (a; = 0, y = 0) . 
The response is measured at several horizontal cross-sections 
below (y = —3, —9) and above (y=6) the source. 



E. Contact force distribution function 

Our simulation allows us to address yet another in- 
teresting and widely-discussed problem: the statistics of 
contact force. Recent experiments indicate that this dis- 
tribution can be well approximated as exponential, that 
is, it is considerably wider than a naively-expected Gaus- 
sian. This is related to the strong heterogeneities of the 
mesoscopic stress in granular matter: it appears to be 
localized to string- like structures known as force chains. 

In the initial sequential packing there is no constraint 
on the sign of the contact force, and its amplitude ap- 
pears to grow indefinitely with the packing depth. After 
the rearrangements, there are no negative forces in the 
system, and therefore their amplitude cannot grow for- 
ever(the total transmitted force is fixed). Figure ^ shows 
the spatial distribution of the contact forces in the sys- 
tems of two different degrees of polydispersity after the 
adaptive stage is completed. One can clearly see that our 
simulations are at least in qualitative agreement with ex- 
periment: it is easy to identify the force chains in both 
cases. 




FIG. 7. Pattern of contact forces from two bead packs of 
different polydispersity. Upper picture has ten percent varia- 
tion in bead radius; lower picture has three hundred percent 
variation. Each contact force is shown as a line joining the 
centers of the two contacting beads. Thickness of the line 
is proportional to the magnitude of the force. In the lower 
picture the bead positions are indicated in light gray. 

We were also able to make a quantitative comparison 
between the simulations and experiments. Figure shows 
the probability distribution function of the contact force 
taken from our simulations of almost monodisperse sys- 
tem. It apparently agrees with the exponential histogram 
observed experimentally. An insight into the origin of 
this exponential behavior is given by the "q-model" due 
to S. Coppersmith et al n^. The further discussion of 
this intriguing result will be published elsewhere |20(| 




FIG. 8. Probability distribution function of the contact 
force in granular packing. 



V. CONCLUSION 

In the study of granular materials, clearcut confirma- 
tion of theories has been elusive. One predicted fea- 
ture of great interest is the null-stress constitutive law 
postulated by Wittmer et a/ ||]. We have verified that 
null-stress behavior occurs in a simplified granular sys- 
tem embodying disorder, perfect rigidity and cohesionless 
contacts. We have measured the free parameter in the 
null-stress law for several situations. We have confirmed 
the validity of our major assumptions used for micro- 
scopic foundation of this constitutive law pq | . Our sim- 
ulation also allowed us to compute directly the ensemble- 
averaged response function, thus providing an additional 
check for the adequacy of the null-stress approach. 

The simulation method has further interesting fea- 
tures. It demonstrates that stable configurations of 



isostatic force networks can be found without chang- 
ing the positions of the nodes. It also reveals order- 
unity variability in the microscopic force distribution re- 
sulting from the relaxation process. Finally, it shows 
strongly heterogeneous response to point forces — much 
stronger than that of the geometric contact network. 
This suggests strong multiple-scattering features in the 
force propagation. 

The observed exponential probability distribution 
function for the contact force is in a good agreement 
with the experiments. Since there is also an indirect ex- 
perimental support for the null-stress law, our choice of 
the system (hard frictionless spheres) appears to be an 
adequate simplification to capture the basic physics of 
granular rigidity. The further simplifications, such as the 
fixed-geometry adaptive algorithm provide an effective 
tool for the future studies of this problem. This may 
include a study of non-linear response of the system to 
large localized perturbation, effects of polydispersity, and 
history dependence of the response. 
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